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We present a precise solution of the polaron problem by 
a novel Monte Carlo method. Basing on conventional dia- 
grammatic expansion for the Green function of the polaron, 
G(k, r), we construct a process of generating continuous ran- 
dom variables k and r, with the distribution function exactly 
coinciding with G(k, r). The polaron spectrum is extracted 
from the asymptotic behavior of the Green function. We com- 
pare our results for the polaron energy with the variational 
treatment of Feynman, and for the first time present precise 
dispersion curve which features an ending point at finite mo- 
mentum. 

PACS numbers: 71.38.+L 02.70.Lq, 05.20.-y 

The polaron problem has a very long history starting 
from the work by Landau |jj . In the most general form it 
is formulated as what happens to the particle when it is 
coupled to the environment, and what are the properties 
of the resulting object, called polaron, which consists of 
the bare particle dressed by environmental excitations. 
This problem arises over and over again not only be- 
cause it is of fundamental importance both for the high- 
energy and for the condensed matter physics, but also 
because the notion of what we call "particles" becomes 
more diverse and new kinds of environments appear. In 
this paper we describe how the polaron problem can be 
solved numerically without systematic errrrs using dia- 
grammatic Monte Carlo method, and present the solution 
of the notorious Frohlich model (see, e.g., Refs. 
First, we explain in detail how this model fits into the 
general Monte Carlo scheme dealing with distribu- 
tion functions of continuous variables. We then describe 
the procedure of extracting the polaron spectrum, E{k), 
from the asymptotic decay of the Green function. Al- 
though for small electron-phonon couplings the polaron 
energy, Eq, and the effective mass, m*, i.e., the bottom 
of the polaron band, are rather well given by the pertur- 
bation theory, the perturbative approach fails to describe 
the spectrum near the threshold E(k) — Eq ~ lo p , where 
Lu p is the frequency of the optical phonon. In fact, the 
threshold features an ending point MM 
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analogous to the ending point of the excitation spectrum 
in 4 He described by Pitaevskii ||]. Our numerical data 
unambiguously confirm this conclusion. 

We start by considering the underlying mathematics. 
Suppose that for a certain random variable/set of vari- 
ables, y, the distribution function ,Q(y), is given in terms 



of a series of integrals with ever increasing number of in- 
tegration variables: 

Q{y) = ^2 X! / dxi ' " dXm F (&n,y,xi, . . .,X m ) . 
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Here £ m indexes different terms of the same order m. 
The term m = is understood as a certain function of y. 
In Refs. |^,|| it was shown how to arrange a Metropolis- 
type stochastic process simulating the distribution Q(y) 
exactly. The process has very much in common with 
the Monte Carlo simulation of a distribution given by a 
multi-dimensional integral. Nevertheless, there is an es- 
sential difference associated with the fact that integration 
multiplicity in the expansion Eq. (|^) is varying. 



FIG. 1. A typical diagram contributing to the polaron 
Green function. 

Projection onto the polaron problem is as follows. Let 
us interpret the Matsubara (imaginary time) Green func- 
tion of the polaron in the momentum-time representa- 
tion, G(k, t), as the distribution function for the random 
variables k and r. We thus identify G with Q, and (k, r) 
with y. Equation (^) is then identified with the dia- 
grammatic expansion of G(k, r) in terms of free-electron 
and phonon propagators within the framework of con- 
ventional Matsubara technique at T = 0. Then, the vari- 
ables Xi,X2, ■ ■ • , x m are the internal times and indepen- 
dent momenta of the diagram £ TO . A typical diagram is 
presented in Fig. 1. Solid lines denote the free-electron 
propagators, G (0) (p, r 2 -t{) = exp[-(p 2 /2-^)(r 2 -n)], 
where /i is the chemical potential. (Plank's constant and 
electron mass are set equal to unity). Dashed lines and 
points stand for phonon propagators, D(q, ti — t%), and 
vertexes of the electron-phonon coupling, V^(q), respec- 
tively. Without loss of generality, one may fix the left end 
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of the diagram at the origin of imaginary time, ascribing 
thus the time r to the right end. 

In this paper we confine ourselves to the Frohlich model 
Pj where phonons are considered to be dispersionless, 
and the electron-phonon coupling has the form 



k,q 



1/2 1 



(3) 
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In Eq. (||) ak and 6 q are the annihilation operators for 
the electron with momentum k and for the phonon with 
momentum q, respectively; a is a dimensionless coupling 
constant. In the Frohlich model the phonon propagator is 
independent of momentum: -D(q, t 2 — 7~i) = exp[— w p (t 2 — 
ti)]. It is convenient, however, to attribute the vertex 
factors to the dashed lines, so that a dashed line with 
the momentum q contributes the factor Z?(q, t 2 — t\) = 
|y(q)| D(t 2 -ti) to the dia gram. The function F is thus 
expressed as a product of G^'s and -D's, in accordance 
with the standard diagrammatic rules. 

Simulating the distribution Q(y) is the process of se- 
quential stochastic generation of diagrams, identical to 
functions F. (In our case these are the diagrams like 
that in Fig. 1, with certain fixed times and momenta.) 
The global process is constituted by a number of ele- 
mentary sub-processes falling into two qualitatively dif- 
ferent classes: (I) those which do not change the type 
of the diagram (change the values of variables of corre- 
sponding function F, but not the function itself), and (II) 
those which do change the structure of the diagram. The 
processes of the class I are rather straightforward, being 
identical to those of simulating continuous distribution 
corresponding to the given function F. In this paper we 
use only one process of this type, namely, shifting in time 
the right end of the diagram Fig. 1. 

In the heart of the method are the sub-processes of 
type II. The generic rules for 
constructing them are as follows Suppose a certain 
sub-process A transforms a diagram F(^ m , y, x\, . . . , x m ) 
into F(£ m+n ,y,xx, . . .,x m ,x m+ i, . . .,x m+n ), and, corre- 
spondingly, its counterpart B performs the inverse trans- 
formation. For n new variables we introduce vector no- 
tation: x — {x m +i, x m +2, ■ ■ ■ , Xm+n}- The process A 
involves two steps. First, it proposes a change, selecting 
a new type of diagram, £ m +„, and a particular value of x. 
The vector x is selected with a certain distribution func- 
tion W(x). There are no requirements strictly fixing the 
form of W(x), but to render the algorithm most efficient, 
it is desirable that W(x) be chosen in accordance with 
some a priori knowledge of coarse-grained statistics of the 
vector x. Upon proposing the modification, the process 
A either accepts it, with probability, P acc (x), or rejects. 
The process B either accepts a modification, removing 



variables x (with a certain probability P lem (x)), or re- 
jects this modification. For the pair of complimentary 
sup-processes to be balanced, the following Metropolis- 
like prescription should be fulfilled 0: 



Race \E) 
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R{x)/W{x), 
1 , 



W{x)/R(x), 
1 , 



if R(x) < W(x) 
otherwise , 



if R(x) > W{x) 
otherwise , 



where 



r>r-*\ PB F {£m+m Ut Xi, ■ ■ ■ , X m , x) 

R(X) = — r 

PA F[^m,y,Xl,...,X m ) 
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and pa and ps are the probabilities of addressing to the 
sub-processes A and B, which, in principle, may differ. 

For solving the polaron problem it is sufficient to have 
only one pair of complementary processes of type II: the 
sub-process A adding a new phonon propagator to the di- 
agram, and its counterpart B removing one phonon prop- 
agator from the diagram. 

Consider the algorithm for the process A. First we 
select the position t\ for the left-hand end of the extra 
phonon propagator. This is done by choosing at ran- 
dom (with equal probabilities) one of the free-electron 
propagators, and by taking for t± any time (with equal 
probability density) within this propagator. Then we se- 
lect the position r 2 for the right-hand end of the phonon 
propagator, in accordance with the distribution function 
cx exp[— u> p (t2— ti)]. After that, we select the momentum 
for this propagator, using the distribution cx (l + q/qo)~ 2 , 



where 



Now the proposing stage is completed, 



and we are ready to perform accept/reject step, follow- 
ing the above prescription, Eq. (^). The corresponding 
function W(x) [x = {ri,r 2 ,q}) reads 



W(x) 



to (1 + q/qo) 2 



,-Up(T2-Tl) 
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where To is the length of the free-electron propagator, 
where the point t% is selected. As mentioned earlier, this 
form of W is by no means the unique one. Apart from 
the factor pg jpj± which will be discussed later, the ratio 
(0) is now completely defined. 

Consider now the algorithm for the process B. We 
simply select at random (with equal probabilities) some 
phonon propagator and with the probabilities given in 
Eqs. (§), (§) remove it. 

To complete the description of the sub-processes A and 
B, we should define the ratio pb/pa- If is quite rea- 
sonable to address to the creation and annihilation pro- 
cedures with equal probabilities. At the first glance it 
might seem that this immediately leads to pb/pa = T 
but this is not true. The point is that when we select an 
electron propagator for placing the point n , we have N e 



2 



equal chances, where N e is the number of free-electron 
propagators in the diagram being modified [denominator 
of Eq. (0) ], and when we select a phonon propagator 
for removing, we have N p h equal chances, where N p h is 
the number of phonon propagators in the diagram from 
which we try to remove the propagator [numerator of Eq. 
(fjj) ] . These N e and N p h are straightforwardly related to 
each other: 

N ph = (N e + l)/2 . (9) 

We thus get 

PB = N e + 1 = N ph 

p A 2N e 2N ph - 1 • [ ' 

As for the processes of type I, these may include 
(i) selection of the time r anywhere on the interval 
(t2n h >°°) according to the simple exponential distribu- 
tion of G^ ' (k, T — T2N ph ) [obviously, the role of the chem- 
ical potential is to make this distribution normalizable, 
and the whole diagram not to diverge to r — > oo; in fact, 
we use /i as a tuning parameter to probe different time- 
scales since the typical length of the diagram in time is 
controlled by the inverse of E(k) — fi ], and (ii) the change 
of the diagram momentum from k to k + p according to 
the distribution function exp[— (k + p) 2 r/2m], where k 
is the average electron momentum of the diagram, i.e., 
k = t -1 J Q dr' k(r'). We find it more convenient how- 
ever, to select the incoming momentum at will and keep 
it fixed, since in this case we collect all the statistics to 
the value of k we are interested in, instead of spreading 
it over the entire fc-histogram. 




Imaginary time 

FIG. 2. Polaron Green function G(k — 0, r) for a — 2 and 
/i = —2.2. Solid line is the exponential fit. 

In Fig. 2 we show the typical data for the polaron 
Green function. Following Ref. ||, we use energy units 



such that u) p = 1. After initial drop at short times we ob- 
serve a pure exponential decay of G(k, r) at longer times 
(provided we are below the threshold of Cherenkov radia- 
tion, E(k) — E < uj p , so that the polaron state is stable). 
From the exponential asymptotic of the Green function 
we readily extract the polaron energy: 

G{k, t > uj- 1 ) Z k exp[-(E(k) - h)t] . (11) 

By fine-tuning the chemical potential very close to E(k) 
we may extend the time-scale for G(k, r) which is given 
by l/(E(k) — /i). Typically, we had reliable statistics 
on the time-scale of order 100/w p , and were thus able 
to deduce the polaron energy to accuracy better than 
O.OlWp. Apart from the polaron energy, the asymptotic 
behavior of the Green function (|ll]) gives us one more im- 
portant physical characteristic of the polaron, the factor 
Zk , which shows the fraction of the bare-electron state in 
the true eigenstate of the polaron: 

Zk = | (free particle fe | polaron k )\ 2 . (12) 




-7.0 



-8.0 

0.0 1.0 _ 2.0 ,. 3.0 4.0 , 5.0 6.0 

Coupling strength 

FIG. 3. Polaron energy Eq as a function of coupling 
strength. Solid line is the Feynman's variational result. 

In Fig. 3 we present our results for the bottom of the 
band Eq as a function of the coupling strength a in the 
most interesting intermediate region < a < 6. As ex- 
pected, our precise data are below the solid line which 
gives the upper bound for Eo (known to be the lowest 
ever obtained for this problem) as derived from the Feyn- 
man's variational treatment || . We cannot but note the 
remarkable accuracy of the Feynman's approach to the 
polaron energy. 

However, the most interesting and instructive data are 
for the polaron spectrum at relatively large k. The per- 
turbation theory result for the dispersion law, E(k) « 
k 2 /2 - a(\/2/ife)sin _1 (A;/\/2) (the solid line in Fig. 4), 
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clearly demonstrates that the first-order correction is sin- 
gular near the threshold of emitting the optical phonon 
and even develops an unphysical maximum (by assump- 
tion, the threshold point was defined as E{k c ) = Eq+uj p ; 
the maximum on the dispersion curve at k < k c is in con- 
tradiction with this assumption) . One is bound to admit 
then that near the threshold the perturbation theory for 
the Frohlich model fails at any a because of the singu- 
lar phonon density of states, which is ^-functional when 
one ignores the curvature of the phonon dispersion law 
u>p(q) ~ ujp = const j^J^] The formalism dealing with 
such cases was developed by Pitaevskii || for the ending 
point in 4 He (a similar approach based on the Tamm- 
Dankoff approximation was suggested in Refs. [^B and 
developed further in HQ). By applying it to the Frohlich 
model we arrive at the following equation for the disper- 
sion law 



u) — a(k — k c ) + bCu 



dx 



x 



+ R(k-k c ,u>) = , (13) 



where Co = w — (Eq+lo p ), R is a smooth function of k — k c 
and Co, and a and b are some coefficients depending on 
a and k c . This equation features an ending point at 
k c , with the parabolic dependence Eq. (Eh at k < k c - 
The Monte Carlo data obtained for a = 1 are shown in 
Fig. 4. We see how an almost perfect agreement with the 
perturbation theory for the band bottom transforms into 
non-perturbative behavior near k c predicted by Eq. (EJ). 



E(k) 




FIG. 4. Polaron dispersion law for a = 1. Solid line is the 
first-order perturbation theory result. 

Apparently, the ending-point is an artifact of the dis- 
persionless phonon spectrum. With the non-zero curva- 
ture of Lo p (q) being taken into account, the ending point 
will transform into sharp crossover from zero to finite 
damping of the polaron state. We estimate the crossover 
region as Ak/k c ~ m/M where M is the mass of the 
host-lattice atoms. 



In summary, we have presented the solution of the 
polaron problem by the diagrammatic quantum Monte 
Carlo method. The method directly simulates the po- 
laron Green function (in the 4D momentum-time contin- 
uum), dealing with its conventional diagrammatic expan- 
sion. The polaron spectrum, extracted from the asymp- 
totic behavior of the Green function, demonstrates es- 
sentially non-perturbative behavior at sufficiently large 
momenta. 

It is worth noting that diagrammatic Monte Carlo ap- 
proach applies to any model dealing with one/few degrees 
of freedom, either continuous or discrete, coupled to the 
thermal bath. Series of the form Eq. (|J) naturally rep- 
resent the partition function in the interaction picture 
f|,||, and this approach was used recently to calculate 
the smearing of the Coulomb staircase in quantum dots 
at strong tunneling conductance jlO) . More generally, di- 
agrammatic Monte Carlo solves any problem which can 
be reduced to Eq. (||). The efficiency, however, severely 
depends on the sign problem, and the convergence be- 
comes very poor if F-functions are not positive definite. 
It is thus of crucial importance to work in the represen- 
tation in which the sign problem is absent. 
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